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ABSTRACT 

We calculate rates of Type I migration of protoplanets in a non-isothermal three-dimensional protoplanetary 
disk, building upon planet-disk models developed in previous work. We find that including the vertical thick- 
ness of the disk results in a decrease in the Type I migration rate by a factor of ~ 2 from a two-dimensional 
disk. The vertical temperature variation has only a modest effect on migration rates since the torques at the 
midplane are weighted heavily both because the density and the perturbing potential are maximized at the mid- 
plane. However, temperature perturbations resulting from shadowing and illumination at the disk's surface can 
decrease the migration rate by up to another factor of 2 for planets at the gap-opening threshold at distances 
where viscous heating is minimal. This would help to resolve the timescale mismatch between the standard 
core-accretion scenario for planet formation and the survival of planets, and could help explain some of the 
rich diversity of planetary systems already observed. 

Subject headings: planetary systems: formation — planetary systems: protoplanetary disks - accretion disks 



1. INTRODUCTION 

The issue of planetary migration due to tidal interactions 
between t he planet and a circum stellar disk was first ad- 
dressed by Goldreich & Tremaine ( 1980), who found that the 
potential of an orbiting planet excites Lindblad resonances 
in the disk, resulting in the transfer of angular momentum. 
They found that the resulting torques on the planet are simi- 
lar in magnitude, so that whether the planet moves inward or 
outward from angular momentum exchange depends o n how 
the torques from the inner and outer disk balance out. IWardl 
(1986) found that given a decreasing radial temperature pro- 
file, the torques from the outer disk are greater than those from 
the inner disk, resulting in a net torque that forces the planet 
inward. Planetary migration via this mechanism is referred to 
as Type I migration to distinguish it from Type II migration, 
which is the radial drift of planets massive enough to clear 
all material out of the Lindblad resonances, thereby forming 
a gap in the disk and moving toge ther with the disk as it is 
accreted onto the star (Ward 1997g|). 

Subsequent work on two-dimensional disks confirms the 
prediction that planets in protoplanetary disks migrate inward 
on short timesc ales on the order of < 10 5 yr for Earth mass 
plane t s (e.g.. lArtymowic z| |1993d: iKorycanskv & Pollackl 
119931 IWardl Il997al) . iTanaka etaiT ( 120021 hereafter TTW) 
find that including the three-dimensional density structure 
of an isothermal disk increases the migration timescale 
by a factor of 2 or 3 ove r a two-dimensional calculation. 
Menou & Goodman (2004) calculate migration rates in a disk 
with a radial temperature profile expected in T Tauri disks 
and also find that the migration timescales increase when the 
thickness of the disk is accounted for. In a disk model that 
includes MHD turbulence, the density fluctuations near the 
planet can be sufficient to change the magnitude and even the 
sign of the total net torque on a small planet, so that the planet 
migrates as a "random walk" rather than a steady inward drift 
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( iNelson & Papaloizoul2004l) . 

While the effect of the density structure in the vicinity of 
a protoplanet on its migration rate has been studied in de- 
tail both analytically and in numerical simulati ons, the ef- 
fect of local temperature variation s has not (e.g., iKlevll 199*91 
iKlev et all2 001: Bat e et al.l2003l) . Most simulations of plan- 
ets embedded in disks employ a vertically isothermal equa- 
tion of state, so that the disk has a simple radially dependent 
temperature structure. Similarly, the analytic calculations dis- 
cussed above assume vertical isothermality for simplicity. In 
this paper, we analyze Type I migration torques from Lind- 
blad resonances in a protoplanetary disk where viscous heat- 
ing and stellar illumination are the primary heating sources. 
Further, we examine how temperature variations induced by 
the presence of the planet itself affect Type I migration. 

Disk temperature variations can affect the net Lindblad 
torques by changing the local pressure gradient in the disk. 
A pressure gradient shifts the angular frequency of the gas 
orbiting the star with respect to Keplerian velocities, which 
in turn shifts the positions of the Lindblad resonances ("Warl 
1997a). The forcing function depends on the distance to 
the planet, so the torques on the planet vary with the loca- 
tions of the Lindblad resonances. To date, no one has ad- 
dressed how a non-isothermal temperature structure will af- 
fect the magnitude of the angular momentum transferred at 
the Lindblad resonances and how that might affect migration 
rates. A variati on in temperature structure , such as those de- 
scri bed in[j ang-Condell & Sasselov)( l2003l hereafter Paper I) 
and Jang-Condell & Sasselov (2004, hereafter Paper II), due 
to radiative transfer effects, may affect the balance of migra- 
tion torques. 

In this paper, we address how the temperature structure of 
the disk affects the locations of Lindblad resonances and the 
consequences for Type I migration. In ^2]we derive the equa- 
tions used for calculating the torques on the planet from Lind- 
blad resonances in a non-isothermal disk. In 50 we describe 
our adopted model disk and briefly discuss how the presence 
of a planet changes the temperature structure of the disk. We 
present our results in |4]and discuss their implications in |5] 
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Finally, in ^6]we present our conclusions. 

2. TORQUES AT LINDBLAD RESONANCES 

In a two-dimensional disk, only one wave mode is excited 
at each Lindblad resonance. In a three-dimensional disk, 
vertical velocity perturbations are allowed, so an entire se- 
ries of wave modes is excited at each Lindblad resonanc e 
(TTW, iTakeuchi & Mivamalfi998t iLubow & QgTIvTe1 1l998). 
The lowest wave m ode, referred to as t he n = mode in TTW 
and the /-mode in Lubow & Ogilvie ( 1998), has no vertical 
velocity component in a vertically isothermal disk, so it is of- 
ten referred to as the two-dimensional wave mode. This cor- 
responds to the wave that is excited in the two-dimensional 
disk, although its properties are somewhat different since it 
is propagated through a disk with a finite height instead of 
an infinitesimally thin disk. The bulk of the angular momen- 
tum, even in therm ally stratified disks, is c arried by the two- 
dimensional wave ( Lubow & Qgilviell998l) . so we shall con- 
centrate our analysis on this wave mode. Vertical motions do 
occur in a thermally stratified disk, but we assume that these 
are small compared to the horizontal velocity perturbations. 
As it is the total angular momentum that the wave carries that 
interests us, the exact form of the wave propagation is outside 
the scope of this paper. 

The pressure gradient is important in the calculation of rel- 
ative torques because it determines the location of the reso- 
nances, indirectly contributing to the magnitude of the torque. 
Local temperature variations can shift the locations of the 
Lindblad resonances by changing the local pressure gradient. 
This so-called pressure buffer can significantly affect the total 
net torques from Lindblad resonances ( Wardl [l997al) . We as- 
sume that the density p and temperature T are not necessarily 
constant with either radius or height, therefore the pressure p 
and sound speed c vary locally as well. 

2.1. Resonance Locations 
Euler's equation for the unperturbed disk is 
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where we adopt an ideal gas equation of state so that p = pc 2 = 
pkT I pm H where k is the Boltzmann constant, p is the mean 
molecular weight (p = 2 for molecular hydrogen), and m H is 
an atomic mass unit. We adopt a cylindrical coordinate sys- 
tem, r = (r, 9,z) where r is the distance from the z-axis and 9 is 
the azimuthal angle. Assuming that the gas in the disk travels 
in circular orbits with v = ril9, then the radial component of 
equation becomes 
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where G is the gravitational constant, and is the mass of 
the star. The equation for the angular velocity of the disk is 
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where fix = y/GM*/r. Defining k = -d\ogp/d\ogr, I = 
-dlog p / dlogT , and h = c/£Ik, we can rewrite equation (0 

as 
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omitting terms of 0(z 4 /r 4 ) and higher. This equation differs 
from calculations of the angular velocity in two-dimensional 
disks which use integrated quantities in that O varies with 
disk height (e.g.. lWardll99 7a). If k a = —d log a/d log r where 
a = j_ pdz is the surface density, then in general, k^k a . 
In the example of a vertically isothermal disk where p = 
J^ h exp(-z 2 /2/i 2 ), the density gradient is 



which yields 
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as in TTW. The angular velocity is constant with z only in 
an isothermal disk, where / = 0. Realistic protoplanetary disk 
models, however, are not isothermal. 

We now consider a planet orbiting the star in the plane of 
the disk at r = a and with angular velocity il p = y/GM^/a. 
The effective location of the mth order Lindblad resonance is 
where r satisfies 
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where k is the local epicycle frequency of the disk, 
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(Artvmowicz l99.3bi lWardl19"88i rr997ai. We assume that k w 
fl, since the departure from a Keplerian flow is small, and 
we find that using equation ^ to solve exactly for n yields 
very similar results for the resonance positions and torques. 
Defining £ = mc/rn, then 



where 



+ 1 , for outer resonances 
— 1, for inner resonances. 
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Then we can solve for the location of the Lindblad resonance, 
a r = r/a, as 



c?J 2 = 



(ii) 



We so lve fo r a r and £ iteratively using equations @ and ( II Q , 
as in iWardi ( 1199 7a). The iteration is necessary because the 
magnitude of the torque is sensitive to the precise position of 
the resonance. Simply calculating a r to first order in h/r can 
introduce a systematic shift in the resonance positions which 
will affect the overall balance of torques. 

2.2. Disk Torque Strength 

In a two-dimensional disk, the magnitude of the torque ex- 
erted near the mth order Lindblad resonance is 



■K 2 ma^ 2 m 
rdD^/dr 



(12) 



where a is the unperturbed surface density of the disk. The 
forcing function of the planet, \l/ m , is defined by Ward ( 1997a) 
to be 

_ rdc/) m /dr+2mf(j) m 

*m — ; ^ j 



TYPE I MIGRATION IN A NON-ISOTHERMAL PROTOPLANETARY DISK 



3 



where / = m(il- f2 p )/f2, and <p m is the amplitude of the mth 
order Fourier component of the disturbing function, 
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where fr™ 2 (a r ) is known as the Laplace coefficient, defined as 
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At the effective location of the Lindblad resonance, / = 

To truly calculate the torques at Lindblad resonances 
in a three-dimensional disk requires a solution to the lin- 
earized fluid equations in three dimensions. For a poly- 
tropic or isothermal equation of state, the eigenmodes can be 
solved for directly" these cases have been studied by others 
tLubow & Ogilviel [T998. TTW). The temperature structure 
of the disk calculated in our model is neither polytropic nor 
isothermal, so we require a numerical solution. To simplify 
the problem, we assume that the vertical velocity perturba- 
tions are small compared to the horizontal ones so that we can 
treat the disk as the sum of independent layers and add up the 
total torque from these layers, as suggested bv lWardl fl988) 
and Artvmowicz ( 1993b). This relies on the assumption that 
relatively little angular momentum is carried away by vertical 
velocity perturbations, so our calculations should be treated as 
a useful estimate of the three-dimensional Lindblad resonant 
torque rather than a rigorous solution. 

For a disk with vertical extent, the forcing function must 
be modified to account for the vertical averaging of the planet 
potential's disturbing function. We can express the torque on 
the disk at height z over an extent dz as 
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where p, r and D» are evaluated at the location of the Lind- 
blad resonance at height z given the temperature structure at 
that height, and ^> m is calculated with a modified Laplace co- 
efficient 
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where ( = z/r. 

As noted in lGoldreich & Tremamel (ll980). for |a r -l| < 1 
and m ^> 1, most of the contribution to the integral comes 
from 9 <C 1. In other words, local conditions around the 
planet have the most effect on angular momentum transfer 
at Lindblad resonances. This is because most of the angu- 
lar momentum transfer occurs from close encounters between 
gas streamlines and the planet. We can then approximate the 
Laplace coefficient by replacing cos 9 with l-8 2 /2 and ex- 
tending the upper limit of the integral to infinity. Then, 
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Substituting in the Bessel function formulation of by 2 (a r , C)> 
2a rV / (a r -l) 2 +C 2 K2m 



lOLr 

(22) 

When £ = 0, this equatio n reduces to the forcing function cal- 
culated by Ward ( 1997a) for a two-dimensional disk. 

The total torque integrated over the vertical extent of the 
disk is then 



T = 

1 tit — 



rdD^/dr 



P^,„dz. 



(23) 



If we assume that the transferred angular momentum is en- 
tirely dissipated in the disk, then the net torque exerted by the 
planet on the disk for the mth order inner and outer Lindblad 
resonances is 



F =T° 



and the total torque is 



m> 1 



(24) 



(25) 



When r is positive (negative), the planet loses (gains) angular 
momentum, and the planet migrates inward (outward). 

2.3. Migration Rates 

Since we consider planets at or below the gap-opening 
threshold, we assume that the torques exerted on the disk by 
the planet are entirely dissipated in the disk and that the torque 
of the disk on the planet is equal to the torques from the planet 
on the disk. We do not take into account the changes to the 
disk made by the action of the tidal torques. We expect that in 
general, the effect of the torques will be to decrease the den- 
sity of the nearby disk by pushing material out of the Lind- 
blad resonances so our torque estimates are upper limits. In 
the absence of feedback form the disk, the rate of change in 
the planet's angular momentum is equal to the total torque, 
and the rate at which the planet moves is 



; 2 TTy/oTrJa y/m 2 [(l-a r ) 2 + ( 2 ]/a r + 9 2 ^yfoT r 
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where K, is the modified Bessel function of order i and 



A = 
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Negative values of da /dt indicate inward migration, positive 
values indicate outward migration. The migration timescale 
is 



(19) 
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3. DISK MODEL 

The model disk is described in detail in Paper II, both with 
and without planets. In this section we will summarize the 
features of the model disk and the temperature changes caused 
by the presence of planets. When calculating the positions 
and torque strengths of the Lindblad resonances, we do not 
take into account the response of the disk to torques from the 
planets. In particular, we do not address density perturbations 
in the disk such as spiral density waves or gap-clearing for 
the larger planets in our models. These density modifications 
will affect both the locations of the resonances as well as the 
magnitude of the torques, but that calculation is outside the 
scope of this paper. 

3.1. The Unperturbed Disk 

To calculate the unperturbed disk structu re , we adopt 
the formalism developed by fCalvetetal. (1991J) and 
iD'Alessio et"allJl998i fT999). with some simplifying assump- 
tions. The model we employ is one of a gaseous disk rotat- 
ing around a hot young star where dust is the primary opac- 
ity source, and the gas and dust are well mixed. The star is 
no longer actively accreting and has no outflow, so the disk 
is quiescent. The star has temperature T+ = 4000 K, radius 
= 2Rq, and mass M* = 0.5M Q . The disk is accreting at 
M = 10~ 8 MQyr _1 and has a Shakura-Sunyaev v iscosity pa- 
rameter of ass = 0.01 (Shakura & Sunvaev 1973). 

To find the temperature structure of the model disk, we 
calculate radiative transfer through the disk with two heat 
sources: viscous heating and stellar irradiation. Figure^sum- 
marizes the resulting vertical temperature and density profiles. 
The upper and lower panels of each plot shows the tempera- 
ture and density, respectively, of the disk versus scale height 
at 0.5, 1, 2, and 4 AU. The dotted line in the lower pan- 
els show the density profile for a vertically isothermal disk 
with the same midplane temperature and surface density for 
comparison. The vertical line in each plot indicates the loca- 
tion of the photosphere, where the optical depth of the disk 
is 2/3. Viscous heating is dissipated at the midplane, so the 
temperature decreases with increasing z. Stellar irradiation at 
the surface creates a temperature inversion, where the upper- 
most layers are hotter than lower layers. The photosphere is 
approximately the location of the vertical temperature mini- 
mum. 

The viscous temperature goes as r" 3 / 4 while the stellar ra- 
diation temperature goes as r" 1 / 2 , so that close to the star, vis- 
cous heating dominates. This causes the midplane tempera- 
ture to be much higher than even the surface temperature at 
small r. Conversely, at larger radii the midplane temperature 
becomes close to the photosphere temperature. 

3.2. Disk with a Planet - Spatial Temperature Variations 

When we insert a planet into the model disk, its gravita- 
tional potential induces a compression of the disk material 
near it, resulting in a decrement in the density at the disk's 
surface. Thus, an isodensity contour at the height of the pho- 
tosphere takes on the shape of a well. When this well is il- 
luminated by stellar irradiation at grazing incidence, it results 
in cooling in a shadowed region and heating in an exposed 
region. 

The temperature variations caused by shadowing and illu- 
mination effects can change the local pressure gradient, which 
will in turn shift the locations of the Lindblad resonances. 
As noted in 32.21 the integral for the perturbing potential is 



weighted toward # <§C 1 for | ct r — 1 1 <C 1 and m ^> 1, implying 
that disk properties immediately around the planet have the 
most effect on angular momentum transfer at Lindblad res- 
onances. Although the equations for calculating torques at 
Lindblad resonances assume an axisymmetric disk structure, 
we can set quantities such as density and temperature equal to 
the local values rather than azimuthally averaged quantities. 

In order to quantify the effect of temperature variations 
on Lindblad torques, we extract the temperature structure 
through a slice of the planet-disk models calculated in Paper 
II, perpendicular to the planet's orbit and through the position 
of the planet and use this to determine the perturbed pressure 
gradient near the planet. This yields new values for the tidal 
torques on the disk, and correspondingly different migration 
rates. Our calculations are meant to put limits on the effect 
temperature variations from radiative transfer effects might 
have on migration rates of planets embedded in disks. 

We use a spline interpolation to calculate the temperature 
between the gridded points, since we want both T and dT / dr 
to be continuous and smooth. Since the density structure 
within about rn of the planet is not well described because 
this material is being accreted onto the planet, the tempera- 
ture structure in this region is also not well understood in our 
models. For this reason, we assume that the temperature sim- 
ply varies linearly across the region within ru from the planet. 
This assumption does not significantly affect the calculation 
of the tidal torques, since the effective resonance locations are 
at \r-a\>2h/3. 

4. RESULTS 

In Paper II, we calculated the three-dimensional temper- 
ature structure around planets embedded in protoplanetary 
disks, taking shadowing and illumination effects at the surface 
of the disk into account. We shall put limits on the effect that 
these temperature variations can have on migration by taking 
a radial slice through the calculated volume at the planet's po- 
sition and finding the migration rates associated with these 
temperature profiles. We calculate the torques at the Lindblad 
resonances and integrate in z from the midplane to the pho- 
tosphere. Note that the photosphere is typically about twice 
the scale height of the disk above the midplane. Above the 
photosphere, the disk is very rarefied, more than an order of 
magnitude less in density, so contributions to the total torque 
are not significant. 

Figure summarizes the total net torques calculated for 
varying planet masses at several distances: crosses for 0.5 AU, 
triangles for 1 AU, squares for 2 AU, and stars for 4 AU. The 
torque is normalized to 

r ° = (S) 2f7 °A^ (28) 

Note that this quantity incorporates the scaling with both stel- 
lar and planetary masses. The solid lines show the torques 
for the unperturbed model disk, and the dotted lines indicate 
the torques for a vertically isothermal disk with the same in- 
tegrated surface density with temperature equal to that at the 
midplane. The dashed lines show how the total net torque 
varies with planet mass when we include the temperature vari- 
ations from radiative transfer. 

In order to better understand how temperature perturbations 
alter migration torques, we discuss how the torque on the disk 
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FIG. 1 . — Vertical temperature and density profiles of the unperturbed disk at (a) 0.5, (b) 1, (c) 2, and (d) 4 AU. The vertical line in each plot shows the location 
of the photosphere. The dotted line shows the density profile of a vertically isothermal disk at the temperature of the midplane with the same surface density for 
comparison. 



varies with z. Noting that 
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dT" 



V dz 



dz 



dz 



we define the scaled torque density as 



, . cto dT tot (to ^ dT° m ' 
p dz p ' dz 

r r m>l 



dT 1 ' 



dz 



(29) 



(30) 



Since 7(2) is scaled by cr , this quantity has the units of torque. 

We plot the scaled torque density in terms of Tq versus z 
at 1, 2, and 4 AU in the upper panels of Figures [5] |4]and 
13 respectively. In each plot, the vertical line indicates the 



scale height of the disk, ho = co/a£l p where Co is the sound 
speed at the midplane of the disk. The short lines on the left 
edge of the plot show the total integrated torque in terms of 
To for the unperturbed model disk (solid line) and a vertically 
isothermal model [dashed line). The line labeled "2D" shows 
total net torque predicted for a two-dimensional disk with the 
same surface density profile as the model disk, and tempera- 
ture given by the midplane of the disk. Note that this value is 
almost exactly equal to the value of the scaled torque density 
at the midplane for a disk with no planet. The solid lines indi- 
cate the torque density profiles for the disk with temperature 
perturbations induced by planets of varying mass, as labeled 
on the plot. The topmost line, labelled M®, shows the scaled 
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FIG. 2. — Total net torques at 0.5 (crosses), 1 (triangles), 2 (squares) and 
4 (stars) AU. The solid lines indicate the total net torque for the unperturbed 
model disk and the dotted lines indicate the total net torque for an isothermal 
disk at the midplane temperature. The dashed lines indicate the total net 
torque for the model disk with a planet of the mass indicated by the horizontal 
axis. 
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FIG. 3. — Scaled torque density profiles, and density and temperature gradi- 
ents as a function of height at 1 AU. The vertical line indicates the disk's scale 
height. (Upper panel): The scaled torque density, in terms of To, is plotted 
versus height from the midplane for a disk with and without a planet, the 
mass of the planet indicated by the labels. The dashed line indicates the net 
torque in a vertically isothermal disk at the midplane temperature. The short 
horizontal lines on left edge of the plot show the total integrated torque for 
the unperturbed (solid line) and vertically isothermal (dashed line) disk mod- 
els. (Lower panel): The temperature gradient, / is plotted as a dashed line. 
The quantities [k + ?,(zSl K / c) 2 /2] (solid line) and [k a + 3/2 + (/ 1 , =0 )z 2 /2ftg] 
(dotted line) are also plotted. 



torque density for the unperturbed model disk. For compar- 
ison, we have plotted the scaled torque density profile for a 
vertically isothermal disk with T = T(z = 0) as a dashed line. 

Close to the midplane, the torque density falls off rapidly 
with z. It then levels off and rises again toward high altitudes. 
The sharp drop-off occurs because of the dependence on z of 
the perturbing potential, indicating that the thickness of the 
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FIG. 4. — Scaled torque density profiles, and density and temperature gra- 
dients as a function of height at 2 AU. See Figure |3]for a full description of 
this plot. 
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— Scaled torque density profiles, and density and temperature gra- 
a function of height at 4 AU. See Figure|3]for a full description of 



disk does play an important role in diluting angular momen- 
tum transfer at Lindblad resonances. All the profiles follow 
the same basic shape near the midplane, including the verti- 
cally isothermal one. 

The torque density for the model disk increases at high al- 
titudes above that of the isothermal disk because of several 
attributes of the disk that all increase the relative differential 
torque higher i n the disk. Firstly, the Lindblad to rque scales 
roughly as c" 2 (Ward 1986; Korvcanskv & Pollack 1993) and 
the vertical temperature minimum roughly corresponds to the 
photosphere, which is the limit of our integration in z. Sec- 
ondly, the radial temperature gradient, I, rises slightly be- 
tween the disk scale height and the photosphere. We plot / 
versus z as a dashed line in the lower panel in each of Fig- 
ures lUHl and |5] The higher values of I favor torques from 
the outer Lindblad resonances by shifting their locations, so 
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planet mass (M e ) 

FIG. 6. — Comparison between migration timescales in the disk-planet 
system in the presence of temperature perturbations and the absence of tem- 
perature perturbations. We plot the migration timescale versus planet mass 
at 0.5 (dotted line), 1 (dashed line), 2 (dot-dashed line) and 4 (solid line) 
AU, showing the perturbed and unperturbed times as lines with and without 
points. (The six-pointed stars represent calculated models, the crosses are 
interpolations for illustration.) 




1 

distance (AU) 

FIG. 7. — Migration rates for the model disk-planet system scaled by the 
planet mass. The upper solid line labelled "2D" is the migratio n rates f or an 
equivalent two-dimensional disk using the formulation of Ward 1 1997a). The 
lower lines are the migration rates calculated for the full three-dimensional 
disk including temperature perturbations from a planet with mass indicated in 
the legend. The solid line, marked Mq in the legend, shows migration rates 
for the disk in the absence of any planet-induced temperature perturbations. 



the net torque rises (IWardlll997al) . Thirdly, the density gra- 
dient, k, is greater than the isothermal density gradient. As 
shown in equation (0}, the location of the resonance is shifted 
by a factor proportional to (3z 2 /2h 2 + k + 1). The quantity 
[fc+3(zJljf/c) 2 /2] in the model disk is plotted as a solid line 
in the lower panel of Figures [3] |4] and |5] For comparison, 
[k a + 3 /2 + (l\ z= o)z 2 /2h 2 ) ], which is the equivalent quantity in a 
vertically isothermal disk, is also plotted as a dotted line. The 
solid line rises above the dotted line most noticeably at 1 AU 
in Figure|5] which is also where the rise in the torque density 
is most significant. 

The total integrated net torques at 1 and 2 AU in the model 
disk are greater than the total torques for an equivalent isother- 
mal disk. This is because the temperature decreases with z up 
to the photosphere in the model disk, so the density falls off 
more than in a vertically isothermal disk, as shown in Figure 
This leads to greater sampling of the torque at the disk 
midplane relative to the photosphere. As the distance from 
the star increases, the vertical temperature profile flattens, and 
the density profile rises relative to the profile of the isothermal 
disk. Contributions to the total torque from the upper layers 
rise relative to contributions from the midplane, effectively 
softening the total torque as a whole, and the total net torque 
of the model disk falls below that of the isothermal disk. 

At 1 AU and inwards, the effect of planet-induced tem- 
perature variations are small. However, as the distance in- 
creases to 2 and 4 AU, planet-induced temperature variations 
become more effective at changing the torque from Lindblad 
resonances. These results are shown in Figure [6] where we 
plot migration time scales versus planet mass at 0.5, 1, 2, and 
4 AU as dotted, dashed, dot-dashed, and solid lines, respec- 
tively. The lines with and without points show migration time 
scales in the model disk with and without considering local 
temperature perturbations, respectively. At 0.5 AU, the mi- 
gration timescales are changed very little by the temperature 
perturbations. Further from the star, the departures from the 



unperturbed disk become more and more pronounced. This 
is, in part, due to the vertical density profile, which weights 
the torques at the midplane more heavily at closer distances 
to the star, as discussed above. But as can be seen from Fig- 
ures |3] |4] and |3 the torques at the midplane are also more 
affected with increasing distance. This is because the total op- 
tical depth of the disk decreases with distance as the surface 
density falls off. This increases shadowing and illumination 
effects at the midplane, which has a greater effect on the net 
Lindblad torques there. 

We re-plot these results in terms of migration rate per Earth 
mass versus distance from the star in Figure0 Given the grid 
of planet distances and masses we have sampled, we interpo- 
late to find the migration rate for selected planet masses and 
plot the variation of migration rate with distance for different 
planet masses, as indicated by the different lines. In the ab- 
sence of perturbations to the disk, the migration rate scales 
with the mass of the planet, so this is good way to concisely 
summarize the effects of local temperature variations. Close 
to the star, even planets at the gap-opening threshold reduce 
the migration rate by only a few percent. However, as the 
distance increases, planet-induced temperature variations be- 
come more and more effective at lowering the migration rate, 
up to 43% at 4 AU. For comparison, we also plot the mi- 
grations rates for a comparable two-dimensional disk, which 
gives migration rates of about twice as fast as the unperturbed 
three-dimensional disk. For this model disk, the overall trend 
is that planetary migration slows as the planet moves closer to 
the star rather than staying constant, although this decelera- 
tion is not sufficient to keep planets from migrating right into 
the star. 
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5. DISCUSSION 

Although the density of the disk is highest at the midplane, 
the contributions to the Lindblad torques at the midplane do 
not solely determine the total net migration torque. The finite 
thickness of the disk is also important because it both modifies 
the perturbing potential and shifts the resonance positions. 

For our model disk, we find that including the thickness of 
the disk reduces the migration rate by about a factor of two 
from a two-dimensional disk. By comparing the model disk 
with an equivalent isothermal disk, we find that including the 
vertical softening of the perturbing potential is, in fact, the 
primary effect, and that the vertical temperature structure has 
a modest effect. We show how the torques in these model 
disks compare to each other in Figure [3] The model disk is 
shown as a solid line, the two-dimensional disk is the short- 
dashed line, and the isothermal disk as a dotted line. 
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FIG. 8. — Comparison of migration torques between different disk mod- 
els vs. radius. All disk models have the same midplane density and surface 
density as the model disk. We show the three-dimensional model disk (solid 
line), a vertically isothermal model (dotted line), and a two-dimensional disk 
(short dashed line). We also show for comparison the calculation done by 
TTW for a three-dimensional (long dashed line) and two-dimensional (dot- 
dashed line) disk. 

TTW also find that the net Lindblad torques in their cal- 
culation of a three-dimensional isothermal disk are a factor 
of ~ 2 - 3 less that of the equivalent two-dimensional disk. 
These torque calculations are also plotted in Figure[8]as a long 
dashed line for the three-dimensional disk and a dot-dashed 
line for the two-dimensional disk. While there is a mismatch 
between the TTW models and our model disk, the results are 
still qualitatively the same in showing how including the ver- 
tical structure of the disk affects the total torque on the planet. 

Planet-induced temperature variations can reduce the mi- 
gration rate relative to that of the unperturbed disk by up to a 
factor of about 2 for planets at the gap-opening threshold at 
around 4 AU. Migration rates are less affected closer to the 
star because of increased viscous heating at the midplane. In 
general, the midplane of the disk contributes most to the res- 
onant torques since temperature variations are minimal both 
because of the high optical depth to the stellar irradiation and 
viscous heating at the midplane. Viscous heating and stellar 
irradiation behave as r" 3 / 4 and r" 1 / 2 , respectively, so viscous 
heating increases faster with decreasing radius than stellar ir- 



radiation. Midplane temperature perturbations are therefore 
washed out by viscous heating closer to the star. This is why 
the presence of a planet has little effect at 0.5 AU, but is highly 

effective at 4 AU. 

The results of Menou & Goodman (2004) also show that 
the thickness of a three-dimensional disk accounts for a sig- 
nificant increase in migration timescales over those of two- 
dimensional disks. They consider an a-disk model similar to 
ours, and find that opacity changes in the disk create local 
regions of increased migration timescale, or decreased migra- 
tion rate, because of local variations in the temperature gra- 
dient. The spikes in the migration timescale that result from 
local changes in the temperature gradient close to the star (< 1 
AU) may be softened by shadowing effects, since these local 
effects may dominate over the large-scale temperature gradi- 
ent. 

As a planet migrates inward, the gap-opening threshold 
mass decreases, so the planet may begin to open a gap and 
undergo Type II migration rather than Type I migration. This 
is a feature of any disk model with a flared structure, where h 
increases faster than r. The gap-opening threshold is where 
rn = (m p /3M ir ) 1 / 3 r = h. Rewriting this equation, the gap- 
opening threshold mass is 

m p = 3M ir f^\ . (31) 

If dlogh/dlogr > 1, then the gap-opening mass increases 
with r. In Type II migration, the planet remains fixed with 
respect to the disk instead of moving through it, as in Type 
I migration. As the disk is accreted onto the star, the planet 
moves with it in the viscous time scale. This is an altogether 
different mechanism for slowing migration. However it still 
suffers from the problem that as long as the disk exists, the 
planet will migrate together with the disk and could still fall 
into the star. 

6. SUMMARY AND CONCLUSIONS 

We have analyzed the effect of local disk temperature vari- 
ations on Type I migration rates of planets embedded in pro- 
toplanetary disks. We include both the vertical structure of 
the disk as well as temperature perturbations caused by the 
presence of the planet itself. 

In comparison to the results of Ward| dl997al) for a two- 
dimensional disk model, the primary difference in migration 
rates is caused by including the vertical thickness of the disk. 
This is because the perturbing potential is diluted by the ver- 
tical distance from the midplane. When we integrate the total 
torque at Lindblad resonances over the vertical extent of the 
disk, we find that this reduces the migration rate by about a 
factor of two. This is consistent with calculations by TTW, 
who find that the vertical dependence of the perturbing po- 
tential reduces migration rates by a factor of ~ 2 - 3 in their 
model disk. 

We estimate the effect of including the vertical temperature 
and density profile of the model disk by comparing it to an 
equivalent vertically isothermal disk and find that the effect 
on the migration rates is modest. This is in part because the 
torque is weighted by disk density, which is highest at the 
midplane, and also because of the rapid decrease in torque 
density with disk height. Therefore disks with the same mid- 
plane temperature and density profiles will yield similar mi- 
gration rates. 

When we include the effects of radiative transfer on the 
perturbation caused by the planet on the disk, we find that 
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the resulting temperature perturbations can further decrease 
the migration rate by up to about a factor of two, depending 
on the planet's mass and position. As the planet increases 
in mass, so does the temperature variation, which results in 
a general decrease in the migration rate. The effect is more 
prominent in planets further from the star, because viscous 
heating decreases faster than stellar illumination with increas- 
ing distance. Viscous heating evens out temperature varia- 
tions, reducing their impact on Lindblad torques. Thus, the 
planet at the gap-opening threshold at the farthest distance in 
our sample is the most effective at reducing its own migration 
rate. 

Our calculation of migration torques given the temperature 
structure in the photosphere of the disk shows a limiting case, 
since we assume that only the local temperature of the disk is 
relevant, and that the azimuthal dependence on temperature is 
unimportant. This is justified by the fact that the perturbing 
potential is strongest closest to the planet, but in principle the 
azimuthal variation of the temperature could dilute the effect 
of local temperature variations on migration rates. 

Another limitation of our planet-disk model is that we have 
not considered the non-linear perturbation to the density struc- 



ture within a few Hill radii of the protoplanet. Since the tidal 
torques from Lindblad resonances scale with disk density, this 
effect may be very important. Numerical simulations univer- 
sally show that spiral structure develops quickly around plan- 
ets embedded in disks. Radiative transfer from stellar illu- 
mination upon these structures can affect the planet and disk 
both directly through heating, and indirectly by changing the 
positions and strengths of the Lindblad resonant torques. 

Still, our model puts interesting constraints on the effect of 
the temperature and density structure on Type I migration. We 
have shown that shadowing and illumination at the surface of 
a protoplanetary disk can slow migration. If the perturbations 
are large enough, they might be able to halt or even reverse 
Type I migratio n, and sav e planet embryos from the "Shiva 
catastrophe" (Ward 1997b). Perhaps this can explain why the 
frequency of giant planets around solar-type stars is relatively 
high, despite the short timescales predicted for Type I migra- 
tion. 



We thank Kristen Menou and our anonymous referee for 
their helpful comments in the preparation of this paper. 
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